Junctions and spiral patterns in generalized Rock-Paper-Scissors models 
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We investigate the population dynamics in generalized Rock-Paper-Scissors models with an ar- 
bitrary number of species N. We show, for the first time, that spiral patterns with iV-arms may 
develop both for odd and even N, in particular in models where a bidirectional predation interaction 
of equal strength between all species is modified to include one N-cyclic predator-prey rule. While 
the former case gives rise to an interface network with Y-type junctions obeying the scaling law 
L oc t 1//2 , where L is the characteristic length of the network and t is the time, the later can lead to 
a population network with iV-armed spiral patterns, having a roughly constant characteristic length 
scale. We explicitly demonstrate the connection between interface junctions and spiral patterns 
in these models and compute the corresponding scaling laws. This work significantly extends the 
results of previous studies of population dynamics and could have profound implications for the 
understanding of biological complexity in systems with a large number of species. 



I. INTRODUCTION 

Non-hierarchical interactions between individuals of 
different species seem to be essential to the development 
of the enormous biodiversity observed in nature. Rock- 
Paper-Scissors (RPS) type models, incorporating some 
of the crucial ingredients associated with the dynamics 
of a network of competing species, are a powerful tool 
in the study of complex biological systems. In its sim- 
plest version, the RPS model describes the evolution of 
3 species which cyclically dominate each other [IH3] (see 
also |4j [5] for the pioneer work by Lotka and Volterra) . 
If the population mobility is small enough, the spatial 
RPS model has been shown to lead to the stable coexis- 
tence of the three species with the formation of complex 
spiralling patterns [TJ [2] . 

The basic interactions behind the RPS game are Mo- 
tion, Reproduction and Predation, but generalizations 
incorporating new interactions and further species have 
also been proposed in the literature [6HT4]. In [11] the 
standard cyclic RPS model was generalized to an arbi- 
trary number of species and it was shown that the emerg- 
ing patterns, in numerical simulations of population dy- 
namics in a cubic grid with periodic boundary conditions, 
depend crucially on whether the total number of species 
is even or odd. While an odd number of species, if the mo- 
bility is not too large, leads to the formation of complex 
spiralling patterns whose characteristic length L remains 
roughly constant in time, an even number of species is as- 
sociated to the formation of partnership domains where 
several species coexist, whose dynamics is controlled by 
interactions of equal strength between species on either 
side of the interfaces separating adjacent domains. This 
has been shown [14] to lead to a scaling regime where the 
characteristic scale of the network L grows as L oc t 1 / 2 , 



which is the typical scaling law associated with grain 
growth and the dynamics of soap froths [T5H22] . 

In the present paper we shall consider a broad fam- 
ily of RPS type models, with an arbitrary number of 
species iV, which may lead to the emergence of spiral 
patterns both for odd and even N. In [14J specific real- 
izations leading to a L ex t 1 ^ 2 scaling law have been inves- 
tigated. This law was shown to accurately describe the 
macroscopic population dynamics of complex networks 
without junctions or with Y type-junctions provided that 
the dynamics is curvature driven, with competing species 
on adjacent interface domains having the same (average) 
strength. In the present paper we shall extend the above 
results to models leading to higher dimensional junctions, 
and investigate the connection between these interface 
junctions and the spiralling patterns which occur under 
cyclic competition between individuals from different do- 
mains. Models with an odd/even number of species shall 
be considered. 



II. GENERALIZED RPS MODELS WITH Z N 
SYMMETRY 

Consider a model where individuals of various species 
and some empty sites are initially distributed randomly 
on a square lattice with AT sites. The different species are 
labelled by i (or j) with z, j = 1, TV, and we make the 
cyclic identification i = i + k N where k is an integer. The 
sum of the number of individuals of the species i (1$) and 
empty sites Ie is equal to J\f. At each time step a random 
individual (active) interacts with one of its four nearest 
neighbors (passive). The unit of time At = 1 is defined 
as the time necessary for M interactions to occur (one 
generation time). The possible interactions are classified 
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Figure 1: (Color online) Scheme of all predator-prey interac- 
tions allowed in a 4-species realization of our model. 

as Motion (active and passive switch their positions) 

Reproduction (active reproduces filling an empty site) 

i — >• ii , 

or Predation (active predates the passive generating an 
empty site) 

i (i — a) —> i or i (i + a) —> i (g , 

where may be any species (i) or an empty site ((g)), 
a = 1, a max with a ma;E AT/2 if N is even or a max = 
(N — l)/2 if TV is odd). 

We denote the corresponding probabilities by mi (Mo- 
tion), n (Reproduction), p^ia (left-handed Predation) 
and pma (right-handed Predation). In this paper we as- 
sume that mi = m, n = r, p Lia = p La and p Ria = p Ra 
for all z, so that all the models have a Zn symmetry. We 
also assume that J\f = 512 2 and that the number densi- 
ties of the various species = Ii/M are all identical at 
the initial time. The number density of empty spaces, 
he = Ie/N, is initially set to be equal to 0.1. Figures 
[T] and [2] show a scheme of all predator-prey interactions 
allowed in 4- and 5-species realizations of our model, re- 
spectively. Except for the labeling of the different species, 
Fig. [I] is invariant under a rotation by an angle of 27r/4 
thus leading to a Z4 symmetry. The same is true for Fig. 
[2] with 27r/4 and Z4 replaced by 2ir/5 and Z5 , respec- 
tively. Full arrows and dashed arrows indicate left- and 
right-handed predations, respectively (with probabilities 
PLa and PRa with a = 1, 2). In the 4-species case i + 2 
or i — 2 coincide and, consequently, Pl2 — PR2 as indi- 
cated by the dashed-dotted arrow in Fig. [I] At each 
time step the active and passive individuals as well as 
the corresponding action are randomly assigned with the 
following probabilities 

i (i + 1) ^\ i (g), i (i- 1) ^ i <8>, 
i (i + 2) ^ i®, i (i-2) ^ i (g, 
i g) — ^ i i, z ^> 

implying that the evolution of a population network is 
stochastic and not fully deterministic. 




Figure 2: (Color online) Scheme of all predator-prey interac- 
tions allowed in a 5-species realization of our model. 



A. Spiral patterns 

Consider a model with N — 4 as a particular ex- 
ample of a model with an even number of species. If 
Pli P and p m = p L2 Pr 2 = (model 14, video 
[23] and top left panel of Fig. [3]), then the partnerships 
{1,3} and {2,4}, between species which do not inter- 
act through Predation, are formed in different spatial 
regions. In contrast, if pl\ = Pl2 = PR2 = P and p R i = 
(model 114, video [24] and top right panel of Fig. [3| 
then there are no viable partnership domains and spi- 
ral patterns with 4 arms do form as a consequence of 
the left-handed predation between neighbouring species. 
For a fixed set of rules, the dynamics of the models may 
be very different for even and odd number of species. 
Consider a model with 5 species (N = 5) with Pli = P 
and pri = Pl2 = PR2 = (model 15, video [25 and 
bottom left panel of Fig. [3|. In this model the species 
tend to organize themselves into spiral patterns whose 
5 arms are dominated by individuals of a species which 
does not interact with individuals of the dominant species 
of the two adjacent arms. For example, individuals of 
the species 1 have as partners individuals of the species 3 
and 4 (note that the colors blue, green, red, yellow, ma- 
genta and white in Figs. [3] and [4] represent the species, 
1, 2, 3, 4, 5 and empty sites, respectively). The fact 
that individuals of the species i predate individuals of 
the species i — 1 (e.g. individuals from species 4 predate 
individuals from species 3) is responsible for the spiral 
patterns in this model, even though predator and prey 
are never the dominant species in adjacent arms. On 
the other hand, if Pli = Pl2 = PR2 — P and pm = 
(model 115, video [26] and bottom right panel of Fig. 
|3| then spiral patterns with 5 arms with a single dom- 
inant species do form, having as immediate neighbours 
individuals from its predator (inside) and prey (outside) 
species. These results can be easily extended to an arbi- 
trary even or odd number of species by considering mod- 
els with p L1 = p and p R1 = p La p Ra (model I) or 
pm = and p L1 p La p Ra p (model II) for a ^ 1. 
The even/odd asymmetry which is present in the case of 
model I completely disappears in the case of model II. 
The results shown in Fig. [3] were obtained with m = 0.5, 




Figure 3: (Color online). Snapshots of simulations of models 
14 (top left), 114 (top right), 15 (bottom left) and 115 (bot- 
tom right). The snapshot of the simulation of model 14 was 
obtained after 1000 generations while the others were taken 
after 20000 generations. 



r = 0.35 and p = 0.15. We have carried out a large 
number of simulations, for a wide range of parameters, 
and verified that the same qualitative results also hold 
for other choices of the parameters m, r and p. 



B. Multiple junctions 

Now let us consider a model with Pli = Pm = p/k 
and Plol— PRa — P for a ^ 1. In this class of models the 
cyclic non-transitive hierarchical predator-prey rule be- 
tween first neighbors does not exist and each species can 
hunt and be chased by any other. In this case N different 
domain types with a single dominant species arise, sep- 
arated by interfaces whose dynamics is curvature driven 
and controlled by interactions of identical strength be- 
tween competing species. Figure [4] shows snapshots of 
4 different simulations with N = 4 and n — 1 (model 

1114, top left and video [27]). N = 4 and k = 20 (model 
IV4, top right and video [28]), TV = 5 and k = 1 (model 

1115, bottom left and video [29]), TV = 5 and n = 20 
(model IV5, bottom right and video [30]). If n = 1 then 
the predation rate Pha — PRa = P does not depend on 
a. In this case, every interface has the same effective 
tension, leading to the formation of an interface network 
with F-type junctions which can be seen in the videos 
[271 [29] and on the left panels of Fig. [3J On the other 
hand, if k > 1, the effective tension of the interfaces sep- 
arating individuals from species i and i db a is not the 
same for a = 1 (smaller effective tension) and a ^ 1 



Figure 4: (Color online). Snapshots of simulations of models 
III4 (top left), IV4 (top right), III5 (bottom left) and IV5 
(bottom right). The snapshots were taken after 5000 genera- 
tions. 



(larger effective tension). This results in the suppression 
of the interfaces with larger effective tensions, and may 
lead to the formation of stable iV-dimensional junctions, 
as shown in videos [28, 30 and on the right panels of 
Fig. [3J As demonstrated in [13] . the average velocity of 
this type of interfaces is proportional to their curvature. 
As a result the network is expected to attain a scaling 
regime where the characteristic length of the network L 
obeys the scaling law L oc t x with A = 1/2. Although 
this has been explicitly demonstrated in [14] , in the case 
of interface networks without junctions or with Y-type 
junctions, we shall provide numerical evidence that the 
same scaling law also applies to interface networks with 
higher-order junctions. All simulations shown in Fig. [3] 
were performed with m = 0.15, r = 0.15 and p — 0.70. 
Although these values were found to be the most ade- 
quate for visualization purposes, we verified that many 
other choices of the parameters m, r and p and k would 
provide similar qualitative results. 



C. Empty spaces distributions 

The black dots in Figs. [5] and [6] represent the distribu- 
tion of the empty sites of the snapshots shown in Figs. 
EJandg] These vacancies are a result of Predation be- 
tween individuals from the competing species. Note that 
in the models 14, III and IV, they are located at the do- 
main borders whereas they are spread out over the entire 
lattice in the case of models 15 and II. Furthermore, the 
number density of empty spaces is much smaller in the 



Figure 5: Spatial distribution of the empty spaces of the snap- 
shots shown in Fig[3] (models 14 (top left), 114 (top right), 15 
(bottom left) and 115 (bottom right)). 



case of model 15 than in the case of model II. This hap- 
pens because in model I the species organize themselves 
to minimize the interactions (other than Motion), while 
in model II adjacent arms are dominated by individuals 
of interacting species. In fact, model II is the simplest 
generalization of the standard RPS model to an arbitrary 
number of species. 



D. Scaling laws 

The average evolution of I^ 1 with time t (over 25 sim- 
ulations with different initial conditions), in the case of 
models 14, III and IV, is shown in Fig. [7| The scaling law 
Ie 1 (x t x describes quite well the late time evolution of 
these population networks, with A = 0.460±0.063 (model 
14), A = 0.479 ± 0.038 (model III4), A = 0.471 ± 0.040 
(model III5), A = 0.483 ± 0.055 (model IV4) and A = 
0.489 ± 0.044 (model IV5). These results were obtained 
considering only the network evolution for t > 100. No 
significant dependence of the scaling exponents on the 
number of species and junction dimensionality have been 
found. In the case of models 14, III and IV empty sites 
appear mainly at the interfaces, whose average thickness 
does not change with time. Hence, the evolution with 
time of the characteristic length L may be estimated as 
L oc I^ 1 . The computed values of A are all very close to 
the A = 1/2, thus extending the validity of the above scal- 
ing law to complex interface networks with junctions of 
dimensionality greater than 3. Note also that the scaling 
regime is achieved earlier in the case of interface networks 



Figure 6: Spatial distribution of the empty spaces of the snap- 
shots shown in Fig. [4] (models III4 (top left), IV4 (top right), 
III5 (bottom left) and IV5 (bottom right)). 

with junctions (models III and IV), a feature also ob- 
served in the case of cosmological domain wall networks 
[31]. The dynamics of spiral patterns in models 15, 114 
and 115 is not curvature driven. After an initial period of 
fast variation the number of empty sites becomes roughly 
constant. This behavior is shown in Fig. [8| where the 
time evolution of I^ 1 is plotted for these models. The 
initial scaling regime is due to the clustering of individu- 
als of the same species to form spiraling domains. After 
some time the characteristic size of the network stabilizes 
as the spiral patterns achieve their stationary size. 



III. CONCLUSIONS 

In this paper we investigated the population dynamics 
in generalized RPS models with an arbitrary number of 
species, demonstrating, for the first time, that: 

1. spiral patterns with N- arms may develop both for 
an odd and even number of species TV; 

2. interface networks with junctions of dimensionality 
greater than 3 may form in models with a sym- 
metric bidirectional predation interaction between 
all the species, with the characteristic length scale 
L of the network obeying the standard scaling law 
L oc t 1 / 2 associated with grain growth and the dy- 
namics of soap froths; 

3. the simplest generalization of the standard RPS 
model to an arbitrary number of species, including 
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Figure 7: (Color online). Time evolution of I E X for models 
14, III and IV. 
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one asymmetric TV-cyclic predator-prey rule, leads 
to a population network containing TV-armed spi- 
ral patterns with a nearly constant characteristic 
length scale at late times. 



This work generalizes the results of earlier investi- 
gations on the evolution of biological populations and 
makes new predictions for the spatial structure and dy- 
namics of population networks whose evolution is de- 
scribed by generalized RPS models. Our work is expected 
to provide a powerful framework for the study of popula- 
tion dynamics and to drive the search for new signatures 
of population dynamics on systems with a large number 
of species. 
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Figure 8: (Color online). Time evolution of i^ 1 for models 15 
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